Genomic copy number variability at the genus, species and population levels impacts in situ ecological analyses of dinoflagellates and harmful algal blooms

The application of meta-barcoding, qPCR, and metagenomics to aquatic eukaryotic microbial communities requires knowledge of genomic copy number variability (CNV). CNV may be particularly relevant to functional genes, impacting dosage and expression, yet little is known of the scale and role of CNV in microbial eukaryotes. Here, we quantify CNV of rRNA and a gene involved in Paralytic Shellfish Toxin (PST) synthesis (sxtA4), in 51 strains of 4 Alexandrium (Dinophyceae) species. Genomes varied up to threefold within species and ~7-fold amongst species, with the largest (A. pacificum, 130 ± 1.3 pg cell−1 /~127 Gbp) in the largest size category of any eukaryote. Genomic copy numbers (GCN) of rRNA varied by 6 orders of magnitude amongst Alexandrium (102– 108 copies cell−1) and were significantly related to genome size. Within the population CNV of rRNA was 2 orders of magnitude (105 – 107 cell−1) in 15 isolates from one population, demonstrating that quantitative data based on rRNA genes needs considerable caution in interpretation, even if validated against locally isolated strains. Despite up to 30 years in laboratory culture, rRNA CNV and genome size variability were not correlated with time in culture. Cell volume was only weakly associated with rRNA GCN (20–22% variance explained across dinoflagellates, 4% in Gonyaulacales). GCN of sxtA4 varied from 0–102 copies cell−1, was significantly related to PSTs (ng cell−1), displaying a gene dosage effect modulating PST production. Our data indicate that in dinoflagellates, a major marine eukaryotic group, low-copy functional genes are more reliable and informative targets for quantification of ecological processes than unstable rRNA genes.

Dinoflagellates encompass most harmful algal bloom (HAB) forming taxa, as well as constituting up to 50% of marine microbial eukaryotic biomass, thus are a major constituent of aquatic microbial ecosystems [13]. Genome size varies considerably in dinoflagellates (~1 Gb to >150 Gb) including some of the largest known eukaryotic genomes, larger than the largest animal (lungfish, 130 Gb) and plant (Paris japonica, 149 Gb) genomes [21][22][23][24][25]. Gene duplication and large-scale expansion appear to have occurred amongst dinoflagellate genomes, and coding genes are often present in multiple tandem repeats [26][27][28][29][30]. Genomes of coral symbiont species (Dinophyceae: Symbiodinaceae) show highly dynamic evolution, driven by gene family expansion via both tandem duplication [28,29] and retroposition [31,32]. Considerable genome size variation and very large genomes occur in multiple planktonic dinoflagellate orders [33], as well as in other groups of marine microbial eukaryotes such as foraminifera, ciliates, and Amoebozoa [33]. GCN of rRNA genes across much of eukaryotic life are considered broadly correlated with genome size [18]. Such large and dynamic genome sizes suggest substantial CNV may exist in these taxa.
Of marine harmful algal blooms forming taxa, those that produce Paralytic Shellfish Toxins (PSTs) are common and have significant public health and economic implications [34]. PST expression can constitute an inducible defence mechanism in marine dinoflagellates in response to the presence of copepod predators [35]. PSTs are synthesised by the cosmopolitan and common marine dinoflagellates Alexandrium species, Pyrodinium bahamense, Gymnodinium catenatum, and Centrodinum punctatum. Dinoflagellate genes associated with PST biosynthesis (sxt) [36][37][38] possess dinoflagellate features such as a unique 22 bp spliced leader sequence on transcripts, a high GC-content, and eukaryotic poly-A tails [36,38]. A relatively low proportion of genes (~10-27%) in dinoflagellates are thought to be regulated at the transcriptional level [21,23,26], with many genes regulated post-transcriptionally. The role of gene dosage acting on this trait may therefore differ in dinoflagellates from that in more highly transcriptionally regulated taxa. Studies of certain species such as A. minutum and A. ostenfeldii have indicated a correlation may exist between cellular PST content and genomic copies of the PST biosynthetic gene sxtA4 [39][40][41]. Some studies have shown that PST synthesis may not be regulated at the transcriptional level [42,43]. This gene has been found to vary in GCN across studies, as A. pacificum and A. catenella show~180-325 copies cell −1 of sxtA4 [40,44,45]; while A. minutum and A. ostenfeldii showed fewer copies, at 1.5-11 copies cell −1 [39,41].The majority of species of Alexandrium had no detectable sxtA4 copies and do not produce PSTs [46]. Thus sxtA4 is a gene with a comparatively lower copy number in dinoflagellates than those that show large scale tandem repeats [26][27][28][29][30]. If consistent across PST-producing species, GCN may constitute a useful marker for in situ ecological analyses of HABs, and potentially for other functional traits governed by genomic dosage.
Because rRNA genes, as compared to coding genes, are likely to be under different selective pressures [47], processes that lead to CNV may differ between them. To determine the impact of CNV on both a functional gene and rRNA barcoding markers, and to examine the role of genome size and time in culture on CNV, we quantified CNV of rRNA genes and sxtA4 in relation to genome size across 51 strains of PST-producing marine dinoflagellate, Alexandrium australiense, A. pacificum, A. catenella and A. minutum. Our selection of strains provided capacity to examine CNV within and between species, in strains maintained in long-term culture, and CNV variance across regions. As diversity analysis employing rRNA genes in particular becomes ubiquitous, we aimed to determine the scale of biases associated with CNV, examine its prevalence across dinoflagellates and indicate potential solutions.

MATERIALS AND METHODS Culture isolation, maintenance and identification
Fifteen non-axenic strains of Alexandrium pacificum were established from a surface net haul collected on 22/11/18 at Mindarie Marina, Western Australia (−31.689127, 115.703103). Single cell isolation of A. pacificum was performed using drawn out glass pipettes and a Nikon Eclipse TS100 inverted microscope (100x magnification). Isolated cells were transferred into Falcon®24 well culture plates containing 1 ml of K/5 medium [48] without sodium silicate. Germanium dioxide was added (5 µg/ml) to prevent diatom growth. Plates were kept at 18°C under a photon flux of 60-100 μmol photons PAR m −2 s −1 with a 12/12 h dark/light cycle (cool white fluorescent). After 3 weeks, the cultures were transferred into 20 ml K media in 70 mL sterile culture flasks (Thermo Fisher Scientific, Massachusetts, USA), and maintained by serial transfer every 3 weeks. In total 36 additional strains of 4 Alexandrium species (Alexandrium catenella, A. minutum, A. pacificum, A. australiense; Supplementary Table 1) were obtained from collections: the Australian National Algae Culture Collection (CS), the Cawthon Institute Culture Collection of Microalgae (CAWD), the Roscoff Culture Collection (RCC), and collections maintained at the Institute for Marine and Antarctic Studies, University of Tasmania. Strains originated from 8 different countries in Europe, Asia, Australasia and the Americas, across states and regions in Australia, and were isolated on differing dates within the past 30 years.
Isolate identity was confirmed by sequencing the D1-D3 region of largesubunit rRNA. Cells were harvested from 50 ml of culture by centrifugation and DNA extracted using the FastDNA spin kit for soil (MP Biomedicals, Santa Ana, California, USA). DNA quality was checked using a Nanodrop 2000 (Thermo Scientific, Waltham, Massachusetts) spectrophotometer. The D1-D3 region of LSU-rRNA was amplified by PCR using primers D1F [49] and D3B [50] in 25 µl reactions containing 5 µl of 5X MyTaq buffer (Bioline, London, UK), MyTaq polymerase (Bioline, London, UK) 0.5 µl, 7.5 pmol of each primer, 1 µg µl −1 BSA (Biolabs, Arundel, Australia), 1 µl of DNA and 15.5 µl of DNA-free water. PCR conditions were 94°C for 5 min, followed by 35 cycles of: 94°C for 30 s, 56°C for 30 s, and 72°C for 1 min; and a final extension step of 3 min. PCR products were verified by 1% agarose gel electrophoresis stained with GelRed (Gene Target Solutions, Dural, Australia) and purified with Zymoclean™ (Zymo Research, California, USA). Sanger sequencing of products was performed by Macrogen (South Korea), with strains assigned to A. pacificum, A. catenella, A. minutum or A. australiense based on comparisons with sequences from verified isolates of each species.

Culture synchronisation and harvest
To measure genome size, CNV and PSTs, Alexandrium spp. were grown to exponential phase in GSe medium [51] at 18°C. For genome size measurement of cell cycle synchronised strains, cells were then incubated for 48 h in darkness to induce synchronisation of cell division [52,53]. Subsamples were fixed with Lugol's iodine at the point of harvest, and cell concentration determined using a Sedgewick-Rafter counting chamber (ProSciTech, Australia) and an inverted light microscope (Leica Microsystems, Wetzlar, Germany). For CNV quantification with qPCR, 60-75 × 10 3 cells in triplicate from each strain were harvested by centrifugation (10 min at 1000 g). For genome size quantification using flow cytometry,~10 5 cells were harvested in triplicate from each strain. At least 10 6 cells were harvested by centrifugation for 10 min at 1000 g for PST measurement.

Genome size measurement
Cells were washed with a 1 × PBS, fixed with 1% (w/v) paraformaldehyde for 10 min, and then washed again with 1 × PBS. Cell pellets were resuspended in 2 mL cold methanol, stored for at least 12 h at 4°C to remove the intracellular chlorophyll, washed twice with 1x PBS, and then stained for > 3 h in 0.1 mg mL −1 propidium iodide and 2 µg mL −1 RNAse (Merck KGaA, Darmstadt, Germany).
A CytoFLEX S Flow cytometer (Beckman Coulter, California, USA) equipped with laser excitation at 488 nm was used for the flow cytometry analysis. BD™ DNA QC Particles Chicken blood cells (3 pg DNA/nuclei; BD Biosciences, San Jose, USA) were used as a standard [54]. In total, 2-µm fluorescent beads were used as stable particles to verify instrument alignment (BD Biosciences, San Jose, USA). Triplicate samples were run at 30 µL min −1 , and data acquired in linear and log modes until at least 10,000 events were measured per sample. Fluorescence emission of propidium iodide stained DNA was detected at 610 ± 10 nm. The peak ratios and coefficient of variation (CV) were quantified with CytExpert software (Beckman Coulter, California, USA). FSC channel was used as a trigger with automatic setting from the manufacturer. Gating and further analysis were only performed for peaks with CV values below 20%, and any peaks above this were rerun. The gating was performed by using FSC-A vs SSC-A gate to exclude debris and use the PI gate on histogram to remove large background noise without the DNA content. Genome size in base pairs used a conversion factor of 1 pg of DNA = 978 Mbp [55].
Genomic copy number quantification with qPCR DNA extraction was carried out using a PowerSoil DNA Extraction kit (QIAGEN, OH, USA) according to manufacturer's instructions. DNA was extracted in triplicate, and quality and quantity determined using a Nanodrop ND-1000 (ThermoFisher Scientific, Waltham, Massachusetts) and R. Ruvindy et al. Qubit 2.0 Fluorometer (ThermoFisher Scientific, Waltham, Massachusetts). qPCR was carried out with a BioRad CFX384 Touch™ System (BioRad, California, USA) using species-specific qPCR assays for Alexandrium rRNA genes [56] and sxtA4 [44] (Supplementary Table 5, Supplementary Fig. 8). qPCRs were run in triplicate using the following cycling parameters: 95°C for 10 s, 35 replicates of 95°C for 15 s and 60°C for 30 s. Total reaction volume was 10 μl, containing 5 μl SybrSelect™ (ThermoFisher Scientific, Massachusetts, USA), 0.5 μM each primer, 1 μl template DNA, and 3 μl PCRgrade water. Samples and Master Mix were loaded to Hard-Shell 384-well PCR plates using epMotion 5075 Liquid Handling Workstations (Eppendorf AG, Hamburg, Germany). Amplification specificity was confirmed using melt-curve analysis. Quantification cycle values were generated by CFX Manager 3.1. Standard curves of sxtA4 and rRNA genes vs quantification cycle (Cq) was developed using ten-fold serial dilution of Gblocks® fragments (Integrated DNA Technologies, Coralville, Iowa, USA). Copy number per µL DNA was determined using the formula: copy number per μl DNA ¼ DNA amount ng ð Þper μl x 6:022 x 10 23 length of fragment 124bp ð Þ x 660 x 10 9 Standard curves, positive controls and negative controls were included in each sample plate as the sample DNA extracted from Alexandrium strains of known concentration (cells µL −1 ). Copies of sxtA4 and rRNA genes per µL −1 DNA were determined relative to qPCR standard curves.

Statistical analyses
The significance of relationships between genome size, rRNA gene copies cell −1 , sxtA4 copies cell −1 and total PST cell −1 were assessed using Spearman's rank correlation and linear regression after transformation, as appropriate, as implemented in GraphPad Prism 7.04. Shapiro-Wilk tests were used to examine normal or log normal distributions. Patterns of genome size and rRNA gene CNV associated with cell culture were based on isolation dates provided by culture collections and isolators. Days from isolation date to sample extraction date were calculated. To account for different laboratory growth rates, cultivation days were converted to estimated number of generations based on published growth rates for each species at culture maintenance conditions (50-80 μmoles PAR; 12:12 L:D cycle18-20°C). Individual strain variance from respective species means were calculated for genome size (pg) and log 10 rRNA gene (copies cell −1 ). Individual strain deviation from species means were calculated using ([(Xi -X̅ )/σ)]. The effect of extended culture periods on variability in genome size and rRNA gene GCN was examined using Levene's test for equal variances among the following three sample groups: (1) Mindarie isolates cultivated for <20 generations; (2) strains cultivated for 100-800 generations; (3) strains cultivated for >1000 generations.

PST quantification
Cell pellets were freeze-dried, then extracted for PSTs and measured [53]. Briefly, samples were resuspended in 2 mL of 1 mM acetic acid and vortexed for 90 s at 100°C for 5 min, sonicated for 5 min and filtered with 0.45 µm PVDF filter (Merck Millipore, Massachusetts, USA). Chromatographic separation (modified [57]) was performed on a Thermo Scientific™ ACCELA™ UPLC system coupled to a Thermo Scientific™ Q Exactive™ (ThermoFisher Scientific, Massachusetts, USA) mass spectrometer. Analysis was performed using an Acquity UPLC BEH
PST content and its relationship to sxtA4 copy number All strains produced PSTs with the exception of A. australiense ATCJ33 and two strains of A. minutum, RC4874 and CCMI1002. The range of PST congeners produced varied considerably between and within species (Fig. 5). Species of A. catenella contained the highest PST content (Figs. 5, 3d) dominated by C1, 2 and GTX1, GTX2, GTX3 and GTX4. Both GTX5, GTX6, as well as decarboxylated versions of STX and NeoSTX were mostly absent. In contrast, A. minutum strains produced a more restricted range of PSTs, with C2, GTX2 and GTX3 produced (Fig. 5). Total PSTs were more uniform among A. pacificum, with no consistently dominant PST variants. A positive relationship was significant between log 10 sxtA4 copies cell −1 and log 10 total PSTs ng cell −1 (Fig. 3d, F = 11.73, df=18, p = 0.003, r 2 = 0.395). A. minutum produced a higher amount of PST per copy of sxtA4 in comparison to A. pacificum. Strains of A. australiense with extremely low or no detectable PSTs (ATCJ33, AADVN1, AT-YC-H), averaged <1 sxtA4 copy cell −1 . In cases where less than one sxtA4 copy cell −1 was found, this was due to copy numbers below the level of quantification of the qPCR assay.

DISCUSSION CNV in microbial eukaryotes and its functional significance
Genomic CNV in eukaryotes, bacteria and archaea has been increasingly documented, sometimes in relation to whole genome duplication or polyploidy [1][2][3][4][5]. CNV is considered a major evolutionary process, influencing the expression of key phenotypic traits [1,2,[4][5][6][7]. The presence of CNV is also of practical significance, impacting the way in which molecular barcoding genes such as rRNA, commonly used for community structure analyses, can be interpreted. The consistency of GCN is poorly understood in microbial eukaryotes, and the causes and consequences of CNV at the genomic level [2,4,5] have been rarely examined. Estimates of CNV in other marine microbial eukaryotes suggested that variation of 1-3 orders of magnitude in rRNA, as was found in foraminifera, ciliates, haptophytes, fungi and diatoms, was considered very high [10][11][12], and was not always related to cell or genome size [10,58]. In this study, rRNA CNV in one genus of a common harmful algae (Alexandrium spp) was found to span~6 orders of magnitude. This should be considered the minimum variation, as we have studied comparatively few species of this genus. Previous estimates of genomic rRNA in Alexandrium spp indicated 10 2 -10 3 copies cell −1 in A. minutum, [16]; from 10 4 -10 6 copies cell −1 in A. pacificum [9], and A. catenella [17,44,45], consistent with our study. Very high rRNA GCN (10 6 ) are also known from other Gonyaulacales spp [59,60] measured using techniques including digital qPCR. Species of Alexandrium therefore span from amongst the lowest to the highest rRNA genomic copies cell −1 of dinoflagellates (Fig. 6a). In our study, most significantly, substantial intra-specific CNV was found among individuals collected from the same population (Fig. 4), of up to 2 orders of magnitude. Recently,~1 order of magnitude CNV among 14 strains of the haptophyte Emiliania huxleyi from different global locations was reported [10]. Using single cell qPCR, two orders of magnitude difference in rRNA gene copies was found within a species of foraminifera [12]. Our results and these recent studies indicate that intraspecific rRNA CNV in marine microbial eukaryotes may be common, suggesting that rRNA gene copies cell −1 of multiple strains is required to determine reasonably representative GCN ranges. While CNV likely impacts analyses of all marine microbial eukaryotes, not addressing this will particularly impact dinoflagellate dominated assemblages, and assuming similarity in GCN in relation to the degree of phylogenetic relatedness appears to be ineffective.
Sources of CNV: genome size variation, time in culture, and phylogenetic divergence Differences in GCN can be the result of genome scale processes including polyploidy, as well as processes impacting part of the genome such as aneuploidy, chromosome duplications, errors during homologous recombination and retroposition [27,31,32]. Mutations and selection while held in laboratory culture could result in CNV, potentially with functional implications, as has been reported in bacteria and fungi, and specific processes may particularly impact rRNA loci [3,61]. In rRNA genes, the presence of tandem repeats that are highly transcribed, and the difficulty of replicating repetitive sequences, make rDNA inherently less stable than other genes and susceptible to CNV [47]. Over cycles of cell division, rRNA GCN in yeast and bacteria in culture showed a pattern of reduction in fungi via recombination-mediated loss [62][63][64]. As well as these inherent genetic factors, environmental factors can lead to induced CNV in rRNA genes in some microbes [58,61,64,65]. In the ciliates Entodinium, Epidinium and Ophryoscolex, GCN of rRNA changes in response to nutrient availability in cultures [65]. Higher rRNA GCN has ecological implications, providing protection against mutagens [62,63] and impacting growth and competitive ability [58,65].
We examined CNV in relation to generations in culture and geographic regions of isolation. Our strains spanned periods in culture of several months to~30 years, yet no significant impact was found on either variability of genome size or rRNA CNV (Fig. 2,  Supplementary Fig. 1), Intraspecific genetic diversity in rRNA gene copies in A. pacificum isolates from sites distant to one another was greater in natural populations than the amount induced by up to 30 years of laboratory selection (Fig. 4). The source of such CNV may be retroposition, as it appears common in dinoflagellates [21,27,31], and was found to be the source of multi-copies of pcna in dinoflagellates [66], as sequences had the remains of dinoflagellate specific spliced leader sequences, which are added to mRNA, evidence of reverse transcription into the genome. In Symbiodiniaceae, chromosomes may be enriched in relation to specific biological processes, suggesting that chromosomes may be duplicated or lost as required, creating CNV [28][29][30]. Whether through a process specific to rRNA, retroposition and/or chromosome duplications, CNV appears to occur over long evolutionary time scales based on our data, as dinoflagellate mutation rates on a decadal time scale were not a significant factor in the high CNV.
It has been calculated that <0.1% of species for which genome size data are available have genomes larger than 100 Gb [24], highlighting the exceptionally large genomes of dinoflagellates [21,22,24,[28][29][30]. In eukaryotes, there is no correlation between genome size and coding sequences in genomes larger than 0.01 Gb, known as the C-value enigma [21,26,67]. Alexandrium species were in the middle to upper half of reported dinoflagellate genome sizes (Fig. 1, Supplementary Fig. 2B), indicating expansion during evolution. Cell cycle synchronised genome sizes of A. pacificum varied~3 fold (Fig. 1), greater than the already large range in previous reports of strains of this species (60-104) pg cell −1 [9,22,40,68,69]. Specialised 'ribosomal chromosomes' were found in the former Alexandrium tamarense species complex, including A. pacificum but not in A. minutum, [68], which in our study and previous studies showed much more consistent genome sizes among isolates (22-29 pg cell −1 in previous studies, as in this study [39,40]). A significant relationship between genome size and rRNA gene copies −1 cell within A. pacificum (Fig. 6a) may indicate that rRNA chromosome duplication contributes to genome size expansion in this species [68], but appears unimportant within the other Alexandrium species examined.
Previous studies have measured dinoflagellate genomes using a variety of methods: flow cytometry, staining with propidium iodide or other DNA dyes, and using cells of chicken or a similarly sized animal genome as a standard, as well as estimates from genome sequencing (i.e. [22, 28-30, 39, 66, 67, 70], Supplementary Table 3). Genome size measurements based on flow cytometry using a standard of an animal genome, which is smaller and differs structurally from a dinoflagellate genome, may lead to larger genome size estimates as compared to estimates of sequenced dinoflagellate genomes (Supplementary Table 3). However, the pattern was not consistent, and some similar estimates were determined regardless of the method (Supplementary Table 3). Dinoflagellates have unusually large genomes, and therefore standards of the appropriate size are generally not available. We ensured that the standard curve used in this measurement Fig. 3 The relationship between the CNV of marker genes and the genome size and total PSTs in Alexandrium species. a CNV of rRNA (±SD) in strains of Alexandrium spp. b Relationship between genome size (pg cell −1 ) and rRNA copies cell −1 in Alexandrium species (F = 22.99, df=49, p < 0.0001, r 2 = 0.319). c CNV of sxtA4 (±SD) in Alexandrium species. d Relationship between total PSTs (ng cell −1 ) and sxtA4 copies cell −1 in Alexandrium species (F = 11.73, df=18, p = 0.003, r 2 = 0.395).
achieved good linearity to minimise any potential error (Supplementary Materials).
Implications for molecular quantification based on rRNA gene barcodes Given the ubiquity of microbial eukaryotic ecological analyses routinely targeting rRNA genes for qPCR and meta-barcoding, it can appear as though such approaches are very well characterised [9,10,13,17,20,56,59,71]. Considerable and justifiable attention has focused on sources of bias that confound universal rRNA gene markers, particularly PCR primer bias, sequencing bias and statistical bias [3,72,73]. Primers have been developed that address primer bias [72,73], statistical approaches have addressed other biases, however less attention has been shown to GCN bias until recently [3,[74][75][76]. CNV is of particular concern for qPCR-based assays that are increasingly applied for the detection of HAB species such as Alexandrium spp [9,16,17,44,45,56,59,60]. Previous estimates of cell abundance using rRNA-targeted qPCR of HAB species indicate both over and under-estimated abundances of up to 500% [45,60,77]. Despite the high variability in rRNA GCN within a population and species (Fig. 4), it is possible to calibrate a qPCR result against another type of environmental cell count via detailed light microscopy or fluorescently labelled flow cytometry [44,45,60]. A sampling design that integrated CNV of multiple populations co-occurring in a region could determine a local rRNA GCN as a calibration proxy. Because mutation rates of rRNA GCN were low over decades (Fig. 2), if a HAB species was seeded via local cyst beds, it may not be necessary to recalibrate with every sample, as doing so would render the simplicity, speed and scalability of qPCR-based quantification obsolete.
Community profiling using molecular barcoding presents a potentially more difficult challenge than qPCR in accommodating CNV, due to unknown variability between taxa and the lack of information on GCN. GCN correction factors are an option for addressing this, and can be incorporated into barcoding analyses [11,[74][75][76]. However, GCN data is very limited and patchily distributed. Current correction approaches rely on prediction using phylogenetically related taxa [65], however GCN of bacterial and fungal rRNA is only moderately phylogenetically conserved [58,65,74] resulting in increasingly inaccurate GCN predictions for diverging lineages [65]. The first study attempting to apply a GCN correction factor for marine microbial rRNA genes used a single mean GCN cell −1 for each class, from a database of~60 species, of 4919 copies cell −1 for dinoflagellates, 166 copies cell −1 for diatoms and 71,710 copies cell −1 for ciliates [77]. While this is an improvement over using GCN as a direct proxy of ASV abundance, the level of variability found in our study shows that such an approach is limited. Genome sizes, potentially a proxy for GCN [ Fig. 6a], are variable and unknown for most dinoflagellates, so this approach also appears to be unfeasible. [65] consider bioinformatic-based solutions to the CNV of rRNA genes is not realistic for microbial eukaryotes, and based on our current data, we agree with that position.
Meta-barcoding of microbial eukaryotes has been referred to as semi-quantitative, as it has been suggested that the relative abundance of taxa may be proportionately representative of the community [20]. However, analyses of bias in metabarcoding of dinoflagellate mock communities using multiple primers targeted to rRNA genes, or comparisons with microscopy counts, found ASV abundances were highly skewed [78,79]. It has been argued that a lack of correlation between cell number and rRNA genes is unimportant, as rRNA gene copies cell −1 is correlated with cell volume [13,19,20,69] in marine eukaryotes, providing a stronger indication of ecosystem function than cell abundance. However, we found only a weak relationship of rRNA gene copies cell −1 with cell volume in dinoflagellates, explaining only 20% of the variance (Fig. 6b, F = 0.26.59, df=106, p < 0.0001, Fig. 4 Copy number variations comparison between A. pacificum isolated from a bloom site against strains from other sites. a Sites where strains of A. pacificum were isolated from the Pacific region including sites across Australia. b Mindarie from where strains of A. pacificum were isolated. c CNV of sxtA4 (±SD) and rRNA (±SD) in Alexandrium pacificum from all sites and from Mindarie only. r 2 = 0.200). This was independent of our present study, as, even when we analysed published rRNA gene copy, cell volume and genome size data [9,10,16,19,22,40,59,60,[67][68][69][70][80][81][82][83][84][85] we found 22% of variance explained (Fig. 6c, F = 16.4, df=58, p = 0.0002, r 2 = 0.2204). When analysing only rRNA copies cell −1 in relation to cell volume for Gonyaulacales spp, the order including many of the most important HAB taxa, we found no significant relationship (Fig. 6c, Supplementary Fig. 2, F = 3.26, df=85, p=ns, r 2 = 0.037). This may be indicative of high evolutionary rates of rRNA genes in this order, noted in relation to long branch lengths in rRNA gene-based phylogenies [86]. Dinoflagellates with the smallest genomes, such as Suessiales spp, appear to show a more straightforward relationship between genome size, cell volume and GCN (Fig. 6a-c), but polyploidy and chromosome duplication are still known from these taxa [85].
Given that ecosystem function is generally of primary interest in microbial ecological research, we suggest metabarcoding be considered a diversity presence/absence measure, and functional genes be quantified using metagenomics or gene specific assays in relation to traits of interest. Our data indicate that comparably fewer copies of the gene related to PST synthesis, sxtA4, were present in PST-producing strains (2-10 2 copies cell −1 ) in Alexandrium, with a smaller range of CNV (Fig. 3c), consistent with studies of sxtA4 in A. minutum (1.5-46 copies cell −1 ; [39,40]) and A. pacificum (34-200 copies cell −1 ; [37,40,44]). Both the reduced CNV, combined with the significant correlation of sxtA4 GCN with PSTs cell −1 indicate that sxtA4 is an informative target for the potential for PST production in situ. We found a variety of PST analogues produced by Alexandrium species, with some species such as A.catenella more consistent in analogues produced than others (Fig. 5). sxtA is involved in the synthesis of the parent compound, saxitoxin, [36,37] and tailoring enzymes in the sxt cluster appear to be responsible for analogues, suggesting that in future these genes could be quantified. An advantage of the sxtA detection is that food web dynamics can be investigated, for example, quantifying sxtA uptake in invertebrates or protists [87]. Examples of similar functional genes that have been or could be detected in situ are cell cycle related genes such as pcna involved in proliferation and growth of HAB species [66], transporters, receptors, genes involved in N and P uptake and other metabolic functions [13,21,29,30,85]. Given the potential for such genes to show a dosage response in dinoflagellates, this could indicate a promising approach to community ecology in dinoflagellate dominated ecosystems.

CONCLUSION
We have shown rRNA gene CNV of up to 6 orders of magnitude in dinoflagellates at the class, genus and species levels, and its relationship to genome size, but not generations in laboratory culture, and very weakly with cell volume. Reported significant correlations between cell volume and rRNA genomic copies in marine microbial eukaryotes were likely driven by taxa with more straightforward genome organisation, rather than dinoflagellates, ciliates and foraminifera, which taken together can constitute the majority of 18 S rRNA signal. We show a significant dosage effect of a functional gene related to a common HAB toxin, and suggest that such low copy functional genes are more stable and informative targets for eukaryotic microbial ecological profiling, using function-based methods such as metagenomics and specific trait-based molecular assays.